#delimit;
clear all;
version 13.1;
pause on;
program drop _all;
capture log close;
set more off;

/*INSERT FOLDER PATH*/
cd /;


set matsize 3500;

use multi_level_dataset_005_field_vars.dta, clear;
drop if ovrlp_anypos==1;
keep id newid pmid year nbcites_it rltd rtrct_year score treat age srce_pmid pubyear yr case_code subv*;
keep if rltd==1;
keep if year==rtrct_year;
xtile qnt_score=score, nq(100);
keep newid qnt_score;
compress;
save interm.dta, replace;

foreach defn in  strct {;

use multi_level_dataset_005_field_vars.dta, clear;
drop if ovrlp_anypos==1;
keep if subv_of_trth_`defn'==1;
keep id newid pmid year nbcites_it rltd rtrct_year score treat age srce_pmid pubyear srce_pmid yr case_code subv*;
merge m:1 newid using interm.dta, keep(match) nogenerate;
compress;

forvalues i = 1(1)100 {;
                gen byte treatXscr`i' = treat*(qnt_score==`i');
                label variable treatXscr`i' "After Retraction � Rltd'ness Score in `i'th pctl.";
                };

label variable treatXscr1 "After Retraction � Rltd'ness Score in 1st pctl.";
label variable treatXscr2 "After Retraction � Rltd'ness Score in 2nd pctl.";
label variable treatXscr3 "After Retraction � Rltd'ness Score in 3rd pctl.";

xtqmlp nbcites_it treatXscr1-treatXscr100 i.age i.yr, fe i(id) cluster(case_code);
estimates save estimates/poisson_main_score_intr_`defn'.ster, replace;

estimates drop _all;
estimates use estimates/poisson_main_score_intr_`defn'.ster;

matrix B=e(b);
matrix V=e(V);

bysort qnt_score: keep if _n==_N;

foreach y of numlist 1/100 {;
	local t=colnumb(B,"treatXscr`y'");
	gen qntl_score`y'_coef = B[1,`t'];
	gen qntl_score`y'_var = V[`t',`t'];
};

	
gen qntl_coefs=.;
gen qntl_vars=.;

foreach j of numlist 1/100 {;
	replace qntl_coefs=qntl_score`j'_coef  if qnt_score==`j';
	replace qntl_vars =qntl_score`j'_var   if qnt_score==`j';
};

gen qntl_95lo=qntl_coefs-1.96*sqrt(qntl_vars);
gen qntl_95hi=qntl_coefs+1.96*sqrt(qntl_vars);

sort qnt_score;

twoway line qntl_95hi qnt_score, lpattern(shortdash) lcolor(red) lwidth(medthin) yline(0, lcolor(black) lwidth(thin)) ||
line qntl_95lo qnt_score, lpattern(shortdash) lcolor(red) lwidth(medthin) ||
line qntl_coefs qnt_score, lpattern(solid) lcolor(blue) lwidth(medium) || if qnt_score >= 1 & qnt_score <= 100, xtitle(" " "100 Percentiles of Relatedness Score") ylabel(-0.75(0.25)0.75, format(%15.2f) angle(horizontal) labsize(small)) xlabel(0(10)100, labsize(small)) legend(off) note(" " "Rank is {it:increasing} in relatedness", span size(vsmall)) saving(graphs/poisson_main_score_intr_`defn'.gph, replace);

};

erase interm.dta;